{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1.6 Error estimation & adaptive refinement\n",
    "\n",
    "\n",
    "In this tutorial, we apply a Zienkiewicz-Zhu type error estimator and run an adaptive loop with these steps:\n",
    "$$\n",
    "\\text{Solve}\\rightarrow\n",
    "\\text{Estimate}\\rightarrow\n",
    "\\text{Mark}\\rightarrow\n",
    "\\text{Refine}\\rightarrow\n",
    "\\text{Solve} \\rightarrow \\ldots\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.geom2d import SplineGeometry\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Geometry\n",
    "\n",
    "The following geometry represents a heated chip embedded in another material that conducts away the heat."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#   point numbers 0, 1, ... 11\n",
    "#   sub-domain numbers (1), (2), (3)\n",
    "#  \n",
    "#\n",
    "#             7-------------6\n",
    "#             |             |\n",
    "#             |     (2)     |\n",
    "#             |             |\n",
    "#      3------4-------------5------2\n",
    "#      |                           |\n",
    "#      |             11            |\n",
    "#      |           /   \\           |\n",
    "#      |         10 (3) 9          |\n",
    "#      |           \\   /     (1)   |\n",
    "#      |             8             |\n",
    "#      |                           |\n",
    "#      0---------------------------1\n",
    "#\n",
    "\n",
    "def MakeGeometry():\n",
    "    geometry = SplineGeometry()\n",
    "    \n",
    "    # point coordinates ...\n",
    "    pnts = [ (0,0), (1,0), (1,0.6), (0,0.6), \\\n",
    "             (0.2,0.6), (0.8,0.6), (0.8,0.8), (0.2,0.8), \\\n",
    "             (0.5,0.15), (0.65,0.3), (0.5,0.45), (0.35,0.3) ]\n",
    "    pnums = [geometry.AppendPoint(*p) for p in pnts]\n",
    "    \n",
    "    # start-point, end-point, boundary-condition, left-domain, right-domain:\n",
    "    lines = [ (0,1,1,1,0), (1,2,2,1,0), (2,5,2,1,0), (5,4,2,1,2), (4,3,2,1,0), (3,0,2,1,0), \\\n",
    "              (5,6,2,2,0), (6,7,2,2,0), (7,4,2,2,0), \\\n",
    "              (8,9,2,3,1), (9,10,2,3,1), (10,11,2,3,1), (11,8,2,3,1) ]\n",
    "        \n",
    "    for p1,p2,bc,left,right in lines:\n",
    "        geometry.Append([\"line\", pnums[p1], pnums[p2]], bc=bc, leftdomain=left, rightdomain=right)\n",
    "\n",
    "    geometry.SetMaterial(1,\"base\")\n",
    "    geometry.SetMaterial(2,\"chip\")\n",
    "    geometry.SetMaterial(3,\"top\")    \n",
    "\n",
    "    return geometry\n",
    "\n",
    "mesh = Mesh(MakeGeometry().GenerateMesh(maxh=0.2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spaces & forms\n",
    "\n",
    "The problem is to find $u$ in $H_{0,D}^1$ satisfying \n",
    "\n",
    "$$\n",
    "\\int_\\Omega \\lambda \\nabla u \\cdot \\nabla v = \\int_\\Omega f v \n",
    "$$\n",
    "\n",
    "for all $v$ in $H_{0,D}^1$. We expect the solution to have singularities due to the nonconvex re-enrant angles and discontinuities in $\\lambda$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = H1(mesh, order=3, dirichlet=[1])\n",
    "u, v = fes.TnT()\n",
    "\n",
    "# one heat conductivity coefficient per sub-domain\n",
    "lam = CoefficientFunction([1, 1000, 10])\n",
    "a = BilinearForm(fes)\n",
    "a += lam*grad(u)*grad(v)*dx\n",
    "\n",
    "# heat-source in inner subdomain\n",
    "f = LinearForm(fes)\n",
    "f += CoefficientFunction([0, 0, 1])*v * dx\n",
    "\n",
    "c = Preconditioner(a, type=\"multigrid\", inverse=\"sparsecholesky\")\n",
    "\n",
    "gfu = GridFunction(fes)\n",
    "Draw (gfu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the linear system is not yet assembled above.\n",
    "\n",
    "### Solve \n",
    "\n",
    "Since we must solve multiple times, we define a function to solve the boundary value problem, where assembly, update, and solve occurs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SolveBVP():\n",
    "    fes.Update()\n",
    "    gfu.Update()\n",
    "    a.Assemble()\n",
    "    f.Assemble()\n",
    "    inv = CGSolver(a.mat, c.mat)\n",
    "    gfu.vec.data = inv * f.vec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SolveBVP()\n",
    "Draw(gfu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimate\n",
    "\n",
    "We implement a gradient-recovery-type error estimator. For this, we need an H(div) space for flux recovery. We must compute the flux  of the computed solution and interpolate it into this H(div) space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "space_flux = HDiv(mesh, order=2)\n",
    "gf_flux = GridFunction(space_flux, \"flux\")\n",
    "\n",
    "flux = lam * grad(gfu)\n",
    "gf_flux.Set(flux)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Element-wise error estimator:** On each element $T$, set \n",
    "\n",
    "$$\n",
    "\\eta_T^2 = \\int_T \\frac{1}{\\lambda} \n",
    "|\\lambda \\nabla u_h - I_h(\\lambda \\nabla u_h) |^2\n",
    "$$\n",
    "\n",
    "where $u_h$ is the computed solution `gfu` and $I_h$ is the interpolation performed by `Set` in NGSolve.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "err = 1/lam*(flux-gf_flux)*(flux-gf_flux)\n",
    "Draw(err, mesh, 'error_representation')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eta2 = Integrate(err, mesh, VOL, element_wise=True)\n",
    "print(eta2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above values, one per element, lead us to identify elements which might have large error.\n",
    "\n",
    "\n",
    "### Mark \n",
    "\n",
    "We mark elements with large error estimator for refinement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "maxerr = max(eta2)\n",
    "print (\"maxerr = \", maxerr)\n",
    "\n",
    "for el in mesh.Elements():\n",
    "    mesh.SetRefinementFlag(el, eta2[el.nr] > 0.25*maxerr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Refine & solve again \n",
    "\n",
    "Refine marked elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh.Refine()\n",
    "SolveBVP()\n",
    "Draw(gfu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Automate the above steps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l = []    # l = list of estimated total error\n",
    "\n",
    "def CalcError():\n",
    "\n",
    "    # compute the flux:\n",
    "    space_flux.Update()      \n",
    "    gf_flux.Update()\n",
    "    flux = lam * grad(gfu)        \n",
    "    gf_flux.Set(flux) \n",
    "    \n",
    "    # compute estimator:\n",
    "    err = 1/lam*(flux-gf_flux)*(flux-gf_flux)\n",
    "    eta2 = Integrate(err, mesh, VOL, element_wise=True)\n",
    "    maxerr = max(eta2)\n",
    "    l.append ((fes.ndof, sqrt(sum(eta2))))\n",
    "    print(\"ndof =\", fes.ndof, \" maxerr =\", maxerr)\n",
    "    \n",
    "    # mark for refinement:\n",
    "    for el in mesh.Elements():\n",
    "        mesh.SetRefinementFlag(el, eta2[el.nr] > 0.25*maxerr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CalcError()\n",
    "mesh.Refine()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the adaptive loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "while fes.ndof < 50000:  \n",
    "    SolveBVP()\n",
    "    Draw(gfu)\n",
    "    CalcError()\n",
    "    mesh.Refine()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot history of adaptive convergence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "plt.xlabel(\"ndof\")\n",
    "plt.ylabel(\"H1 error-estimate\")\n",
    "ndof,err = zip(*l)\n",
    "plt.plot(ndof,err, \"-*\")\n",
    "\n",
    "plt.ion()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
